NASA Contractor Report 3183 


Loan o: •; ■- 

A'r W .L. i . ; 

. ' KiRTLAi-.'L, 



-:r.:r.uR.|j|| 

- 


A Compressible Solution 
of the Navier-Stokes Equations 
for Turbulent Flow About an Airfoil 





S. J. Shamroth and H. J. Gibeling 


CONTRACT NASl-15214 
OCTOBER 1979 




TECH LIBRARY KAFB, NM 


NASA Contractor Report 3183 


TECH LIBRARY KAFB, NM 

1 


□□bis 


A Compressible Solution 

of the Navier-Stokes Equations 

for Turbulent Flow About an Airfoil 


S. J. Shamroth and H. J. Gibeling 
Scientific Research Associates , Inc. 
Glastonbury , Connecticut 


Prepared for 

Langley Research Center 

under Contract NAS 1-1 52 14 


NASA 

National Aeronautics 
and Space Administration 

Scientific and Technical 
Information Branch 


1979 




TABLE OF CONTENTS 


Page 

SUMMARY 1 

INTRODUCTION 2 

LIST OF SYMBOLS 7 

ANALYSIS 10 

The Governing Equations 10 

The Turbulence Model 16 

The Numerical Procedure 22 

The Coordinate System 24 

RESULTS 26 

NACA 0012 Airfoil at Zero Incidence 26 

CONCLUDING REMARKS 34 

APPENDIX - THE NUMERICAL METHOD 36 

Linearization Technique 36 

Application of Alternating-Direction Techniques 41 

Solution of the Implicit Difference Equations 43 

REFERENCES 4 5 

FIGURES 51 


i i i 



SUMMARY 


A compressible time-dependent solution of the Navier-Stokes equations 
including a transition-turbulence model is obtained for the isolated airfoil 
flow field problem. The equations are solved by a consistently split linearized 
block implicit scheme due to Briley and McDonald. A nonorthogonal body fitted 
coordinate system is used which has maximum resolution near the airfoil surface 
and in the region of the airfoil leading edge. The transition-turbulence model 
is based upon the turbulence kinetic energy equation and predicts regions of 
laminar, transitional and turbulent flow. Mean flow field and turbulence field 
results are presented for an NACA 0012 airfoil at zero and nonzero incidence 
angles at Reynolds number up to one million and low subsonic Mach numbers. 



INTRODUCTION 


Flow over an isolated airfoil, particularly in the presence of separation, 
is a problem which clearly demonstrates the importance of developing efficient 
and accurate numerical solutions for the Navier-Stokes equations. The isolated 
airfoil problem arises in a variety of practical applications including the 
aircraft wing, control surfaces, propellers and helicopter rotor blades. When 
these components are at small or modest incidence, accurate flow field 
predictions can be obtained by combining inviscid and boundary layer analyses 
in a noniterative mode (Refs. 1 and 2). If the boundary layer-inviscid flow 
interaction remains localized, accurate flow field predictions can be obtained 
even with small regions of separated flow present (Ref. 3). However, if 
significant regions of separation are present, then the usual boundary layer 
assumption that the pressure distribution about the airfoil is either 
independent of the viscous flow development near the surface or can be 
approximated by a simple inviscid displacement surface, becomes significantly 
in error and a full Navier-Stokes solution which simultaneously predicts the 
pressure and velocity fields is required if valid predictions are to be 
obtained . 

Turning to the helicopter problem in particular, a phenomenon of primary 
interest is that of dynamic stall. Static stall occurs when an airfoil is placed 
at large incidence in a steady stream; dynamic stall occurs when the incidence is 
a function of time. Dynamic stall differs from its static counterpart in two 
major ways. First of all, the maximum lift obtainable under dynamic conditions 
is greater than that under static conditions. Secondly, though under static con- 
ditions lift is uniquely related to incidence, under dynamic conditions stall 
depends upon the time history of motion and has a hysteresis loop associated with 
it. As the helicopter blade travels through the rotor disc, the blade experi- 
ences a varying incidence angle. Over most of the disc the blade will be 
unstalled (i.e., the flow will not contain any large separated regions leading 
to a decrease in blade lift or a generation of large blade moment coefficients); 
however, over a portion of the disc large regions of separated flow may appear 
and over these regions the blade performance will deteriorate. This time history 
dependence of the problem makes dynamic stall prediction a particularly difficult 
task. 
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Despite its complex nature dynamic stall has been the subject of a number 
of theoretical and experimental investigations over the past years. The 
behavior of the leading edge separation bubble has been investigated by both 
Velkhoff, Blaser and Jones (Ref. 4) and Isogai (Ref. 5) and the dynamic stall 
mechanism itself has been scrutinized by McCroskey and Phillipe (Ref. 6), 
McCroskey, Carr and McAlister (Ref. 7) and Parker (Ref. 8). More recent 
investigations have been reported at the AGARD Symposium on Unsteady Aerodynamics 
(Ref. 9). These include the unsteady flow investigation of Saxena, Fejer and 
Morkovin (Ref. 10) and the dynamic stall review paper of Philippe (Ref. 11). 

In addition to the large number of experimental investigations, a variety 
of analytical approaches have focused upon the dynamic stall problem. A review 
of these methods has been given by McCroskey in Ref. 12. Among these approaches 
are inviscid vortex shedding models of Ham and Garelick (Ref. 13) and Bandu, 

Sanger and Souquet (Ref. 14), the data correlation approach of Carta and his 
coworkers (e.g.. Ref. 15) and the semi-empirical models of Ericsson and Reding 
(Ref. 16) and Lang (Ref. 17). These were followed by the approaches of Crimi 
and Reeves (Ref. 18), Shamroth and Kreskovsky (Ref. 2) and Kreskovsky, Shamroth 
and Briley (Ref. 1), all of which are based upon solutions of fluid dynamic 
equations. More recently, motivated by a desire to relieve some of the 
restrictive assumptions and/or empiricism of earlier methods, attention has 
focused upon Navier-Stokes solutions to the isolated airfoil problem. In an 
early work of this type, Mehta and Lavan (Ref. 19) solved a stream function 
vorticity formulation of the laminar incompressible Navier-Stokes equations to 
predict flow about an impulsively started airfoil. Although this method 
required considerable computer run time, its excellent results convincingly 
demonstrated the practical benefits which could be realized from Navier-Stokes 
solutions. In a recent work Mehta (Ref. 20) used a more efficient numerical 
procedure to solve incompressible laminar flow about an airfoil oscillating 
through incidence regimes in which stall occurs. Another Navier-Stokes analysis 
is that of Wu and Sampath (Ref. 21) and Wu , Sampath and Sankar (Ref. 22) who 
have applied the Wu-Thompson integro-dif ferential formulation (Ref. 23) to the 
impulsively started airfoil problem. In a similar vein Kinney and Cielak 
(Refs. 24 and 25) have applied a vorticity formulation of the Navier-Stokes 
equations to the incompressible airfoil flow field problem and Hodge and Stone 
(Ref. 26) have investigated stalled airfoils again using a Navier-Stokes solution. 
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Although arguments can be made In favor of one of these procedures versus 
the other. It is clear that as a group these efforts have demonstrated that 
application of Navier-Stokes formulations to the airfoil problem is both 
feasible and practical. However, these procedures all have had three major 
limitations associated with them: (i) incompressibility, (ii) laminar flow 

and (iii) lack of generality of airfoil geometry. In regard to the first of 
these items, the preceding analyses all are incompressible and none can be 
extended readily to the compressible case. In regard to the second limitation, 
all these analyses assume the flow to be laminar although presumably this 
assumption can be relieved in a straightforward manner if simple eddy viscosity 

and forced transition concepts are accepted. Insofar as geometry is concerned 
it is noted that coordinate systems for the Navier-Stokes solutions are at 
least one degree more stringent in the coordinate system smoothness than those 
required for inviscid analyses, and hence a straightforward extension of an 
existing coordinate system developed for solution of the inviscid equations may 
not be possible. The Navier-Stokes analyses mentioned above have used a 
coordinate system which either is conformal or is generated from a solution of 
Poisson's equation. The conformal requirement is a very stringent one which 
could make it very expensive or even impossible to investigate many airfoil 
shapes of practical interest by virtue of the resulting coordinates having 
numerous and rapid changes in curvature. The requirement that the coordinate 
system be generated from a solution of Poisson's equation is not nearly as 
constraining as is the requirement that it be conformal. In principal most if 
not all airfoils of interest could be generated in this manner. However, as 
discussed by Thompson, Thames and Mastin (Ref. 2 7), in practice convergence 
problems can arise when obtaining a solution to the Poisson equation required 
to generate the coordinate' grid. In addition practical difficulties can arise 
in specifying a coordinate grid having high resolution in regions where the 
greatest need for resolution exists; i.e., in regions in which the dependent 
variables change rapidly. As a result of these restrictions, a more general 
coordinate generation system would be a distinct advantage in applying Navier- 
Stokes analyses to the airfoil flow field. 

The problem of eliminating the incompressible assumption from the full 
Navier-Stokes equations for airfoil flow field calculations has been overcome 
by Verhoff (Ref. 28), Deiwert (Ref. 29), Levy (Ref. 30) and Gibeling, Shamroth 
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and Eiseman (Ref. 31). Verhoff applied MacCormack's fully explicit method 
(Ref. 32) to the airfoil problem; however, since the procedure is fully explicit, 
a small time step is necessary to maintain numerical stability as a result of 
the locally refined mesh in the boundary layer and long computer run times 
result. In this regard conditionally stable schemes such as fully explicit 
schemes are not an optimum choice when mesh refinement is required for boundary 
layer definition; in these schemes the allowable marching step size is limited 
by the spatial step size leading to large run times. On the other hand, 
unconditionally (in a linear sense) stable schemes such as some of the implicit 
schemes do not suffer from this characteristic. Both Deiwert’s (Ref. 29 and 
Levy's (Ref. 30) analyses are based upon MacCormack's more recent hybrid 
implicit-explicit-characteristics scheme (Ref. 33). By virtue of an enlarged 
stability bound this new procedure is much more efficient than the original 
MacCormack procedure (Ref. 32) for airfoil calculations; however, it does 
present formidable coding problems. 

Another approach to the airfoil problem has been taken by Steger (Refs. 

34 and 35) who has used an approximate viscous analysis in conjunction with 
the coordinate generation procedure of Thompson, Thames and Mastin (Ref. 36) 
to predict laminar flow about an airfoil. The equations solved by Steger 
(Refs. 34 and 35) are not the full Navier-Stokes equations but rather a reduced 
set of equations which retain only those viscous terms important in thin shear 
flows. As a consequence Steger's approach could be invalid when the shear 
layer no longer remained thin or even became significantly inclined to the 
coordinate along which the thin shear layer approximation is made. Both of 
these problems could be particularly significant in the stall problem under 
consideration here. The numerical scheme used by Steger is that given by Beam 
and Warming (Ref. 37), who used the time linearization approach of Briley and 
McDonald (Ref. 38) in conjunction with the approximate factorization technique 
of D'Yakonov (see Yanenko, Ref. 39) to produce a linearized block implicit 
scheme. More recently Briley and McDonald (Ref. 40) have shown that the Beam 
and Warming (Ref. 37) derivation produces essentially the same consistently 
split linearized block implicit scheme as that developed earlier by Briley and 
McDonald (Ref. 38) using the Douglas-Gunn ADI procedure (Ref. 41) in its 
natural extension to systems of partial differential equations. Gibeling, 
Shamroth and Eiseman (Ref. 31) applied this Br i 1 ey-McDona Id procedure to solve 
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the Navier-Stokes equations for the flow about a cylinder and an airfoil, 
although in this study, as in Refs. 28, 34, and 35, the calculations were 
performed for laminar flow. 

The present effort extends the compressible flow work of Gibeling, Shamroth 
and Eiseman (Ref. 31) and concentrates upon the latter two of the previously 
mentioned problem areas: (ii) turbulent flow and (iii) coordinate generality. 

The procedure adds a turbulence model to the airfoil flow field and demonstrates 
the practicality of constructive nonorthogonal coordinate systems by producing 
calculations about an NACA 0012 airfoil using a constructive coordinate system. 
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LIST OF SYMBOLS 


Except where dimensions are specified, all quantities in the following 
are nondimens ional ; physical velocities are normalized by u r> density by 
P r , pressure by P r u r , dynamic viscosity by u r » and time by (L/u r ) where L 
is the reference length. 


y 

D e 

ay 

E E 

ay 

F 6 . 

ija 


1J 

ij 


G . . 


Turbulence constant 
Turbulence structural coefficient 
Coefficient in pressure-velocity law 
Coefficient in pressure-velocity law 
Airfoil chord 

Turbulence structural coefficient 
Turbulence structural coefficient 
Turbulence structural coefficient 

Coefficient in tensor form of momentum equation 

Coefficient in tensor form of momentum equation 

Coefficient in tensor form of momentum equation 

Vector in Navier-Stokes equations 
Vector in Navier-Stokes equations 

Metric tensor coefficient 
Inverse metric coefficient 

Vector in Navier-Stokes equations 
Vector in Navier-Stokes equations 

Coefficient in tensor form of momentum equation 
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LIST OF SYMBOLS (CONT’D) 


v 

v 

w 

X 


Coefficient in tensor form of momentum equation 
Jacobian 

Turbulence kinetic energy 
Turbulence length scale 

Coefficient in tensor form of momentum equation 
Reference length, taken as airfoil chord 
Coefficient in tensor form of momentum equation 
Pressure 

Coefficient in tensor form of momentum equation 
Reynolds number 

Coefficient in tensor form of momentum equation 
Turbulence Reynolds number 

Coefficient in tensor form of momentum equation 
Time 

Cartesian velocity 

Reference velocity, usually taken as free stream velocity at 
upstream infinity 

Cartesian velocity 

Contravariant velocity component 

Vector in Navier-Stokes equation 

Cartesian coordinate parallel to chord 
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LIST OF SYMBOLS (CONT'D) 


y 


a 

y 



Cartesian coordinate normal to chord 
Distance normal from surface 
General coordinate 
Boundary layer thickness 

Kronecker delta 


£ 

n 

K 

u 


o 

< 


a 

e 


Turbulence energy dissipation 

- General coordinate 

- Von-Karman Constant 
Viscosity 

Reference viscosity, taken as free stream viscosity at upstream 
inf ini ty 

Turbulent viscosity 
General coordinate 
Density 

Reference density, taken as free stream density at upstream 
inf ini ty 

Prandtl number for turbulence energy transfer 
Prandtl number for turbulence dissipation transfer 


T 


Shear stress 



ANALYSIS 


The Governing Equations 

The Mean Flow Equations 

The flow field about an airfoil is governed by the Navier-Stokes equations 
and in conjunction with a suitable turbulence model a solution of the time- 
dependent form of these equations would serve to predict the flow field for both 
laminar and turbulent flows. The form of the equations expressed in the more 
common coordinate systems can be found in standard fluid dynamic texts (e.g.. 

Ref. 42 and 43) and the equations themselves have been derived in general 
tensor form by McVitte (Ref. 44) for inviscid flow and by Walkden (Ref. 45) for 
viscous flow. 

When solving the Navier-Stokes equations for flow about airfoils, basic 
decisions must be made concerning the coordinate system to be used and the form 
of equations to be solved. The presence of flow field bounding surfaces which 
do not fall on coordinate lines presents significant difficulties in applying 
boundary conditions and unacceptably large truncation errors could result. One 
method of eliminating this problem transforms the governing equations from the 
usual Cartesian coordinates to a new coordinate system having the airfoil 
surface as a coordinate line; in general, the process may lead to a nonorthogonal 
coordinate system. 

Several candidate forms of the governing equations present themselves for 
solution in the general nonorthogonal case. For example, Gibeling, Shamroth 
and Eiseman (Ref. 31) have solved the set of Navier-Stokes equations in general 
tensor form using the density and the contravariant components of thp velocity 
vector as dependent variables. The governing momentum conservation equations 
were the components of the vector equations aligned in the coordinate directions. 
The equations were written in the form 

Continuity Equation 

(^J) + _i_ {pvPj)= 0 (1) 

at 
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Momentum Equation 
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where t is time, J is the Jacobian, g „ is a metric tensor coefficient, v is 

th 

a contravar iant component of velocity in the i coordinate direction, p is 
density, y 1 is a coordinate and 
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and A and B relate pressure, velocity and density through the gas law equation 
which is written for constant total temperature as 

P ■ ^)[a + Bg fj v'v J ] (5) 

ct R 

In the previous equations g are coefficients of the inverse metric tensor, 
y is viscosity. Re is Reynolds number, Latin indices are summed from 1 to 3 
and Greek indices are summed from 0 to 3; v° is defined as unity. In addition 
all quantities have been made nondimensional ; physical velocities have been 

normalized by a reference velocity, u , density by a reference density, P , 

2 r r 

pressure by p^u^. , dynamic viscosity by and time by L/u^ where L is a 

reference length. This represents the first candidate form of the equations. 
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Solution of the equations in this form yielded excellent results for flow 
about a circular cylinder (Ref. 31), however, resolution difficulties resulting 
from large truncation error were encountered when the contravariant velocity 
components were used as dependent variables for flow about an airfoil. Therefore, 
the authors of Ref. 31 solved the airfoil flow field by employing the physical 
velocity components in the coordinate direction as dependent variables rather 
than the contravariant components of the velocity vector. This was accomplished 
by noting that in an orthogonal coordinate system the physical velocity 
component u^ could be expressed in terms of metric data and the contravariant 
velocity components v 1 as. 


Uj = 


( 6 ) 


where no summation on the repeated index is implied. Therefore v 1 was replaced 
in the governing equations by (v 1 / /g” ) /g^ and the equations solved for 

u. = v 1 / / g . . - When this approach was used, a solution for flow about a 

X 11 

Joukowski airfoil was obtained without difficulty, and although Eq. (6) is 
restricted to orthogonal coordinates a more general expression of this type 
could be used for nonor thogonal coordinates. This approach represents the 
second candidate form. 


The third candidate form solves a divergence form of the Navier-Stokes 
equations in which the Cartesian velocity components are the dependent 
variables. The governing equations (for constant total temperature flow) are 
continuity and the momenta equations in the two coordinate directions; however, 
the independent spatial variables are transformed from the Cartesian 'coordinates 
x,y to a new set of coordinates £,n where 


f * €(x,y,t) 

i? « 7?(x ,y,t) 

r-t (7) 
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As shown by Steger (Refs. 34 and 35) the governing equations then become 


dW/D a [ Wf. FF G£ y 1 a I W17. F17 X G-Fjy 

dr d£ [ D D D J di) [ D D D 



where 


( 8 ) 


o • - € y *? x 



Steger used these equations in conjunction with a thin shear layer stress 
approximation to predict an airfoil flow field (Refs. 34 and 35) and this 
approach was carefully considered as a candidate approach in the present 
investigation, although of course modifying the approach so as to retain the 
full stress tensor appropriate to the Navier-Stokes equations. 

The final candidate form considered in the present investigation also is 
based upon density and the Cartesian velocities as dependent variables and 
continuity and the momenta equations written in the Cartesian directions as 
governing equations. However, in this approach the equations are not solved 
in divergence form but rather are solved in the quasi-linear form 


aw dw dF 

JT + JT + 


1 

Re 
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x d£ 


dF, dGj dGj 

+ V *~Frj + ^y“aT + 


dGj 


( 10 ) 
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In deciding upon an approach to use in the present effort, the first 
candidate form was eliminated from consideration because of the resolution 
difficulties it exhibited in the airfoil calculation of Ref. 31. An examination 

O 

of the expressions for I and E p in Eq . (3) shows that the second candidate 

av av 

form requires second derivatives of metric coefficients which in turn requires 
third derivatives of y a with respect to the Cartesian coordinates x and y. 

In contrast the third and fourth candidate approaches only require evaluations 
of second derivatives of y a with respect to the Cartesian coordinates. This 
requirement of one less degree of differentiability upon the grid coordinate 
transformation represents a distinct advantage for the third and fourth 
candidates, since a less smooth coordinate system is required and can be more 
easily generated. For this reason the second approach was eliminated from 
further consideration. Thus the final choice was narrowed to these latter two 
possible approaches. 

Before proceeding to airfoil calculations, both the third and fourth 
procedures were used to repeat the cylinder flow field calculation previously 
presented in Ref. 31. The results of this calculation are shown in Fig. 1. 

The calculations were run until the maximum change in each dependent variable 
during one body traverse time divided by the dependent variable Itself dropped 
to below one per cent. As can be seen in Fig. 1, the results of the fourth 
candidate approach, Eq. (10), were in good agreement with the three sets of 
results presented in Ref. 31. In contrast the results of the third candidate 
approach showed considerable disagreement with all other results. 

A careful examination of the equations was performed to determine the 
reason for the poor predictions obtained from formulation number three. After 
a term by term examination of the equations, the major difficulty was found to 
reside in the momenta equation pressure gradient terms. In this formulation 
the pressure term is of the form 9Gp/3£ where G is a function of the geometric 
data. The numerical representation of this term is 


dGp G<e + AC)p{£ + AC)-G<C-AOp<£-Af) 

d( ' ' 2A? 

( 11 ) 
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In the cylinder calculations, p varies from grid point to grid point in the 
fourth significant figure whereas G varies in the first or second significant 
figure. As a result pressure changes were commensurate with the geometric 
data truncation error, an obviously undesirable situation. 

Although the cylinder represents only one case, the resulting poor 
treatment of the pressure gradient term by the third formulation is still 
expected to be a problem in performing airfoil calculations. The problem 
would be somewhat alleviated with increasing Mach number since increasing Mach 
number is accompanied by increasing pressure difference from point to point. 

In addition, the cylinder geometry represents a relatively stringent test of 
formulation three. The airfoil is a much more streamlined body than the 
cylinder with the coordinate curvature being much less pronounced over most of 
the flow field. Nevertheless, large coordinate curvature does occur in the 
airfoil leading edge region. Therefore even if the problems exhibited by the 
strong conservation approach are not as serious in the airfoil case as they 
proved to be in the cylinder case, they are still expected to be present. A 
more detailed discussion of the problem is given by Shamroth and Gibeling in 
Ref. 46 . 

Based upon these considerations, the fourth coordinate formulation was 
chosen for the present effort. 


The Turbulence Model 

Since the present effort addresses the problem of turbulent flow, it is 
necessary to specify a turbulence model suitable for this problem. One compli- 
cating factor in hypothesising and applying a turbulence model for the isolated 
airfoil flow field is that the flow is not turbulent everywhere. Far from the 
airfoil the flow is laminar. In addition, even near the airfoil surface the 
flow is laminar in the vicinity of the airfoil leading edge. Thus any proposed 
model must be capable of dealing with laminar, transitional and turbulent flow. 
The turbulent flow problem including laminar and transitional regions has been 
treated approximately but quite successfully for the simpler problem of boundary 
layer flows by McDonald and Fish (Ref. 47), Shamroth and McDonald (Ref. 48) and 
Kreskovsky, Shamroth and McDonald (Ref. 49). These investigators applied an 
integral form of the turbulence kinetic energy equation in conjunction with the 
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boundary layer mean flow momentum equations to predict a wide variety of flows 
in both forward and reverse transition. As shown in Refs. 47-49, the model gives 
excellent agreement between prediction and data for a variety of typical test 
cases. In an alternative but similar approach Jones and Launder (Ref. 50) and 
Launder and Jones (Ref. 51) used a two-equation turbulence model in predicting 
relaminarizing flows, but did not apply this model to the forward transition 
problem. Application of the model to the forward transition process in boundary 
layers was initiated by Pridden (Ref. 52) who obtained predictions of transitional 
boundary layers (see Launder and Spalding (Ref. 53)). Later the Launder-Jones 
model was applied to the forward transition problem again in boundary layers 
by Quemard and Michel (Ref. 54); however, none of the calculations presented 
in Ref. 54 proceeded successfully through transition. In summary, at the 
initiation of the present effort, calculations through transition in both the 
forward and reverse direction had been made routinely for boundary layer flows 
using an integral form of the turbulence kinetic energy equation (Refs. 47-49); 
however, questions remained concerning the problem of predicting transitional 
flows via the differential equation of turbulence kinetic energy for airfoil 
flow fields. 

The approach taken in the present effort assumes an isotropic turbulent 
viscosity, y^, relating the Reynolds' stress tensor to mean flow gradients. 




/ au i 

dU: N 
+ — 


du k 


\dx ] 

dXj y 

3 

d* k 

5 'll 


( 12 ) 


Using Favre averaging (Ref. 55) the governing equations then are identical to 
the laminar equations with velocity and density being taken as mean variables 
and viscosity being taken as the sum of the molecular viscosity, y, and the 
turbulent viscosity, y^. Originally, it was intended that the turbulent 
viscosity be obtained from the two-equation turbulence model in which the 
turbulence kinetic energy, k, and turbulence dissipation rate, e, were taken as 
dependent variables. This model has been used by several investigators (e.g.. 
Refs. 50 and 51) for fully turbulent flow fields and, in fact, has been 
utilized successfully by Gibeling, McDonald and Briley (Ref. 56) in a study of 
turbulent combusting flows made using a reacting flow version of the same 
numerical procedure applied in the present investigation. 
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The equations governing the development of k and e have been given by 
several authors (e.g.. Launder and Spalding (Ref. 53 )) and in Cartesian 
coordinates are given as 


dpk dp uk 

dt + dx 


+ 




dpuk d C f fjL T \ dk 

*~ir * 

du k \ d Uj dk dk 

ax^ ) dx k ' 


(13a) 


d/5€ dpU€ d/5V€ 
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( a ' Ui 'll 

^ <* x k 


/ 

C 2 k 2 / X Mt 



(13b) 


In Eqs . (I3a) and (13b) and o £ are Prandtl numbers for k and £ respectively 
and and ere empirical functions. The equations are discussed in some 
detail in Ref. 53 and this discussion will not be repeated here. The turbulent 
viscosity is related to k and e via the Prandtl-Kolmogorof constitutive equation 

= P c n k * /e 


(14) 


In the present analysis the following values were assumed 

cr ( * 13 
cr k * 1.0 
C, = 1.55 


(15) 
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For fully developed turbulent flow Cy = .09 and 0 .^ = 2.; in relaminarizing flows 
Jones and Launder (Ref. 50) give 


= 0-09exp [- 2.5( I. + R t /50. )] 
C 2 - 2-0 {l. 0-0-3 [ exp(-R^)]} 


however, although this expression has given good results in relaminarizing 
flows, it has not led to the successful prediction of flows undergoing forward 
transition in Ref. 54. Therefore, it was felt if a transition model were to be 
obtained, Eq. (16) should be modified and the successful integral turbulence 
energy procedure of (Refs. 47-49) was used for guidance. This procedure 
utilizes a turbulence function where 


a 


l 



(17) 


and a^ 


is taken as a function of a turbulence Reynolds number of the form 


a , * 

f(R r ) ' 

/ 1.0 + 6.66 Cl 

r f(R T ) i 

— i — i 
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where 


a 0 = .0115 

f( R t ) * 100. r t 0 2Z R t < I 
f(R T ) = 68 I R t + 614.3 R t >40 


( 19 ) 


19 


I 


and a cubic curve was fit for values of R^ between 1 and 40. As previously 
discussed. Ref. 47-49 utilized an integral form of the turbulence kinetic energy 
and therefore R^ was defined as an average value. 



( 20 ) 

In the present effort was defined as the local ratio of turbulent to laminar 
viscosity, a^ was evaluated via Eq . (18), related to a^ via Eq. (17) and 
C 2 was evaluated via Eq. (16). 

In the early stages of the present effort the mean flow Navier-Stokes 

equations were solved in conjunction with the partial differential equations 

governing turbulence kinetic energy and turbulence dissipation equations given 

earlier; however, problems arose with the solution. The problem areas were 

those regions in which the flow was expected to be laminar or transitional 

rather than fully turbulent. Regions of this type are found far from the 

airfoil and in the vicinity of the airfoil leading edge. One problem which 

arose concerned evaluation of turbulent viscosity via Eq. (14). In regions 

where k is small, y should be small; however, in terms of the model y depends 
2 T T 

on the ratio of k /e and therefore regions of small k can exist simultaneously 

with regions of large y^ if the predicted value of e is small enough. Since 

regions of small turbulence kinetic energy should coincide with regions of 

small turbulent viscosity by definition, such behavior is unacceptable. 

The problem of ensuring small turbulent viscosities far from the airfoil 
regardless of the values of k and e was solved by modifying the viscosity law 
so that 





,- e -< v 8>/8 



( 21 ) 


for y g > 6 
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The use of Eq. (21) ensured small turbulent viscosities far from the airfoil. 

The second problem, the region near the airfoil leading edge, proved much more 
troublesome. In this region the calculation led to negative values of k and/or 
e and then became unstable. It is possible that this behavior either represents 
a flaw in the model for laminar and transitional flow or is related to starting 
the solution from an unrealistic initial k-e profile. The physical source of 
the problem has not yet been determined, although the numerical problems clearly 
reside in a stiffness of the governing partial differential equation (McDonald, 
Ref. 57). Because of the difficulties encountered, the two-equation model was 
not pursued further in this effort but rather a model combining the turbulence 
kinetic energy equation with a specified length scale was used. 

In this model the length scale was taken as a minimum value of two lengths; 
a wall length and a wake length. The wall length was assumed to be given by a 
conventional wall damped Prandtl's mixing length, viz 



( 22 ) 

with a maximum value of 0.09 6. In Eq . (22) tc is the von Karman constant taken 
as 0.43, y+ is the dimensionless distance from the airfoil surface and 6 is the 
boundary layer thickness. The wake length scale was taken as i = .056 where 6 
is the wake thickness. Rather than update the boundary layer wake thickness 
at each time step, for the present the analysis assumed the boundary layer 
thickness over the first ten per cent chord region to be given by 


8 - k, x ,/z + k 2 

and over the last eighty percent chord to be given by 


( 23 ) 


s - k » 08 

3 ( 24 ) 

where k^, and k^ were chosen at convenient program restart time steps to fit 
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the generated velocity solution. The two 6 formulations were smoothly joined 
in a blending region, 0.10 f_x/c 0.20. Although a step-by-step update of 6 
based upon the calculated velocity or shear stress field would be preferable, 
in the interim the present approach assures a smooth variation of 6 with 
distance along the chord and is considered acceptable for the present demonstra- 
tion. It is anticipated that in subsequent work the local value of the boundary 
layer thickness will be used and the present strategy was used simply to avoid 
possible numerical problems arising from the determination and use of a time 
instantaneous boundary layer thickness. 

When applying the turbulence kinetic energy equation with a specified length 
scale, the turbulent viscosity appearing in the turbulence kinetic energy source 
term, p (3u . /3x +3u,/3x . ) 3u . /3x, , is replaced by Eq . (21) and the dissipation 

X 1 K K. 1 1 K 

rate e is expressed in terms of k and 2, by 



(25) 

The resulting equation then is solved to predict the turbulence kinetic energy 
f ield . 

The Numerical Procedure 

The numerical procedure used to solve the governing equation is a 
consistently split linearized block implicit scheme originally developed by 
Briley and McDonald (Ref. 39) and embodied in a computer code termed. MINT, 
an acronym for Multi-dimensional Implicit Nonlinear Time-dependent. The basic 
algorithm was further developed and applied to both laminar and turbulent duct 
flows by Briley, McDonald and Gibeling (Ref. 58). Subsequently, a three- 
dimensional compressible Navier-Stokes combustor flow analysis employing the 
MINT procedure was developed by Gibeling, McDonald and Briley (Ref. 56) and 
this procedure was then employed by Levy, Shamroth, Gibeling and McDonald 
(Ref. 59) to determine the feasibility for computing a turbulent shock-wave 
boundary layer interaction with an implicit method. More recently Gibeling, 
Shamroth and Eiseman (Ref. 31) have applied the Briley-McDonald scheme to flow 
about a cylinder and an isolated airfoil. Since the scheme has been outlined 
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in detail in several publications available in the open literature, it will not 
be detailed here. Rather only a brief outline of the procedure will be given 
with more detail in the Appendix; the reader interested in further detail is 
referred to Refs. 31, 38 and 58). 

The method can be outlined as follows: the governing equations are replaced 

by an implicit time difference approximation, optionally a backward difference 
or Crank-Nicolson scheme. Terras involving nonlinearities at the implicit time 
level are linearized by Taylor expansion about the solution at the known time 
level, and spatial difference approximations are introduced. The result is 
a system of multidimensional coupled (but linear) difference equations for the 
dependent variables at the unknown or implicit time level. To solve these 
difference equations, the Douglas-Gunn (Ref. 41) procedure for generating 
alternating-direction implicit (ADI) schemes as perturbations of fundamental 
implicit difference schemes is introduced in its natural extension to systems 
of partial differential equations. This technique leads to systems of coupled 
linear difference equations having narrow block-banded matrix structures which 
can be solved efficiently by standard block-elimination methods. 

The method centers around the use of a formal linearization technique 
adapted for the integration of initial-value problems. The linearization 
technique, which requires an implicit solution procedure, permits the solution 
of coupled nonlinear equations in one space dimension (to the requisite degree 
of accuracy) by a one-step noniterative scheme. Since no iteration is required 
to compute the solution for a single time step, and since only moderate effort 
is required for solution of the implicit difference equations, the method is 
computationally efficient; this efficiency is retained for multidimensional 
problems by using ADI techniques. The method is also economical in terms of 
computer storage, in its present form requiring only two time-levels of storage 
for each dependent variable. Furthermore, # the ADI technique reduces multi- 
dimensional problems to sequences of calculations which are one-dimensional in 
the sense that easily-solved narrow block-banded matrices associated with one- 
dimensional rows of grid points are produced. Consequently, only these 
one-dimensional problems require rapid-access storage at any given stage of 
the solution procedure, and the remaining flow variables can be saved on 
auxiliary storage devices if desired. 


23 


The Coordinate System 


The choice and construction of a coordinate system is an important component 
in the successful solution of the isolated airfoil problem. The presence of 
bounding surfaces in a flow field which do not fall upon coordinate lines 
present significant difficulties for numerical techniques which solve the 
Navier-Stokes equations. If a bounding surface (such as the airfoil surface) 
does not fall upon a coordinate line, undue numerical errors may arise in the 
application of boundary conditions and considerable effort may be required to 
reduce these errors to an acceptable level. Although this problem arises in 
both viscous and inviscid flows, it is more severe in viscous flows where 
no-slip conditions on solid walls can combine with boundary condition truncation 
error to produce numerical solutions which are both qualitatively and quanti- 
tatively in error. Thus coordinate systems are sought in which each no-slip 
surface of the specific problem falls on a coordinate line. 

In addition to requiring that no-slip surfaces fall on coordinate lines, 
other considerations constrain the choice of a suitable coordinate system. 

The two most important of these are smoothness of the coordinate system and 
distribution of grid points. Since the governing equations (e.g., see Eq. (10) 
of the present report) contain derivatives of the transf ormation, any trans- 
formation chosen must have a degree of smoothness and continuity which allows 
accurate representation of these derivatives. Secondly, in viscous flow 
problems different regions have different characteristic length scales and 
thus fine grid resolution must be available in some regions of the flow field 
whereas only relatively coarse resolution need be maintained in other regions. 

For example, in calculating turbulent flow about an airfoil at a chord Reynolds 
6 -4 

number of 10 , grid spacings on the order of 10 chord may be required near 
the airfoil surface whereas spacings of the order of the airfoil chord may be 
permissible in the far field. « 

In considering possible airfoil coordinate systems, several candidate 
schemes are available. Among these candidates are conformal coordinate systems 
such as that used by Mehta (Ref. 20), systems based upon solution of Poisson’s 
equation, and a constructive system (Refs. 31 and 60). Although conformal 
systems are excellent choices when they can be used, their very nature restricts 
their generality and hence they were not considered for the present effort. 


24 



A more general method is that of Thompson and his coworkers (e.g., Thompson, 
Thames and Mastin (Ref. 36)) which solves a Poisson equation to generate 
coordinates. Although this approach is considerably more general than the 
conformal approach, the generation of coordinates by this approach requires 
solution of a partial differential equation and obtaining a converged solution 
of Poisson's equation for general geometry can be a problem (Ref. 36). In 
addition difficulties may be encountered in obtaining a desirable distribution 
of grid points. 

The third approach is that of constructing a coordinate system in which 
the required airfoil is by definition a coordinate line and in which grid 
point placement is specified by the user. Such a procedure was developed 
originally for the isolated airfoil problem by Gibeling, Shararoth and Eiseman 
(Ref. 31) and extended to the cascade problem by Eiseman (Ref. 60); it is this 
procedure which was used in the present calculations. To the authors' knowledge 
the present calculations represent the first successful attempt at solving the 
Navier-Stokes equations about an airfoil using a constructive coordinate system. 
The details of the constructive system are presented in Ref. 31 for an isolated 
airfoil and in Ref. 59 for a cascade. 
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RESULTS 


NACA 0012 Airfoil at Zero Incidence 

The first calculation presented is for a NACA 0012 airfoil at zero 
incidence. A highly stretched coordinate system was used for the calculations 
with grid spacing normal to the airfoil being of the order of .0005 c (where 
c is the chord) in the vicinity of the airfoil surface and of the order of 0.8 c 
at the outer grid boundary, approximately some three chords from the airfoil 
surface. The coordinate system was constructed using the method described in 
Refs. 31 and 60. The grid, which consisted of 41 x 30 points, extended 
approximately three chords upstream of the airfoil, three chords above and 
below the airfoil and six chords downstream of the airfoil. The dependent 
variables were the Cartesian velocity components, u and v, the density, p, and 
the turbulence kinetic energy, k. The governing equations were taken as the 
continuity equation, the momenta equations, and the turbulence kinetic energy 
equation. The independent variables for the turbulence kinetic energy equation 
were transformed from the Cartesian coordinates, x and y , in the same manner 
as the momenta and continuity equations were transformed from the Cartesian 
to the computational coordinates. 

The calculation was initiated by specifying the inviscid solution throughout 
the flow field; applying zero velocity conditions on the airfoil surface and 
assuming a turbulence kinetic energy field. Boundary conditions were set as 
follows: on the airfoil surface the Cartesian velocities were set to zero and 

the boundary layer approximate transverse momentum equation (9p/8n = 0) used. 

At the outer boundary, velocity and density were set to the inviscid. values for 
a Joukowski airfoil of approximately the same thickness whereas at the downstream 
boundary, first derivatives were set to zero along the coordinate direction 
approximately aligned with the free stream. A sketch of the boundaries is 
presented in Fig. 2. The turbulence kinetic energy equation was treated by 
setting, k, equal to a free stream turbulence level of 0.001 at the upstream 
outer boundary and setting the (nearly streamwise) derivative of k equal to 
zero at the downstream boundary. At the airfoil surface the correct condition 
on the turbulence kinetic energy should be k=0, however, the turbulence kinetic 
energy reaches a maximum very near the wall and then declines to zero (Ref. 61). 
Rather than try to include enough grid points in the near wall region to define 
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this sharp maximum and the subsequent decrease to zero at the wall an artificial 
boundary condition was applied in the spirit of the so-called wall function 
formulation. This boundary condition takes the form of setting the turbulence 
kinetic energy derivative at the wall equal to zero and thus not defining the 
decrease from the maximum k to the zero k wall condition. Obviously, an 
improved formulation would include enough points in the near wall region to 
adequately define the problem there and such an approach is expected to be 
followed in subsequent efforts. 

The calculation was run in three steps. Initially a spatially varying 
turbulent viscosity field was specified and the mean flow equations were solved 
with these equations. During this part of the calculation since the transients 
were not of interest the equations were solved with the matrix conditioning 
approach developed by McDonald and Briley (Ref. 62) which results in a spatially 
varying time step. Matrix conditioning considerably accelerates the convergence 
of the mean flow equations to steady state. After the major flow readjustment 
had occurred the mean flow was "frozen" and the turbulence equations were solved. 
Finally all equations were solved and the flow development to steady state 
computed. Using this technique the solution reached convergence in approxi- 
mately 150 time steps. Although other factors may be involved this result 
should be contrasted to the 800 time steps required by Steger (Ref. 34) to reach 
convergence in solving the thin shear layer equations. 

A sketch of the predicted surface pressure distribution is presented in Fig. 

3, which also shows the predictions of Mehta (Ref. 20) for laminar incompressible 

flow at Re = 10,000 and the data of Gregory and O’Reilly (Ref. 63) for turbulent 

flow at a Reynolds number of 2.88x10 . In the present calculations the reference 

length, velocity, density, viscosity and temperature were chosen to obtain a 

£ 

Reynolds number based on chord of 10 and a Mach number of approximately 0.2. As 
can be seen in Fig. 3, both predictions are in good agreement with data over the 
forward portion of the airfoil. However, over the last 40% span, the measurements 
indicate substantially more recompression than either analysis predicts and, in 
fact, both analyses show good agreement with each other. 

The suction surface pressure discrepancy between the analyses and the data 
deserves further comment. The trailing edge pressure predictions of Mehta’s 
analysis were made with a coordinate grid which concentrated the available grid 
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points in the trailing edge as well as the leading edge regions. Thus the 
coordinate system used in the present analysis which obtains relatively high 
resolution in the leading edge region at the expense of some resolution in the 
trailing edge region does not seem to be the source of the discrepancy. However, 
it must be recalled that Mehta (Ref. 20) considered a laminar calculation 
whereas the Gregory and O’Reilly (Ref. 63) data is for turbulent flow. A laminar 
boundary layer is more susceptible to separation than a turbulent one and indeed, 
the prediction of Mehta (Ref. 20) does show separation upstream of the trailing 
edge whereas the data show the boundary layer to remain attached over the entire 
suction surface. Therefore the discrepancy between the data of Ref. 63 and the 
analysis of Mehta could be the result of the computed laminar boundary layer 
separating and modifying the trailing edge pressure distribution. 

The present analysis however considers turbulent flow, but due to grid 
point limitations the boundary layer may not be sufficiently resolved for 
optimum results. Thus either the existing boundary layer resolution or perhaps 
the assumed turbulence model could cause the predicted boundary layer displace- 
ment effects to be larger than actually occur and this would result in a lower 
than expected pressure in the suction surface trailing edge region. This 
possible explanation is given some weight as a result of the second calculation 
presented in the present report. In this second calculation the zero incidence 
airfoil flow field was predicted using a grid having somewhat better resolution 
in the boundary layer region and the result was a higher prediction of trailing 
edge suction surface pressure. This will be discussed subsequently. 

Isovel contours in the immediate vicinity of the airfoil surface for the 

zero incidence NACA 0012 case are presented in Fig. 4. The horizontal axis of 

the plot is x/c which is the dimensionless streamwise distance along the airfoil 

centerline. The vertical axis is the normal distance from the airfoil surface 

y /c. A sketch of these quantities is shown in the bottom portion of the 
s 

figure. The contours in the vicinity of the leading edge stagnation point, 
x/c = 0, represent changes in the essentially inviscid velocity along the 
stagnation streamline, however, downstream of the leading edge region, x/c >0.1, 
the isovels represent the viscous flow region and here the growth of the boundary 
layer is evident. The airfoil trailing edge is located at x/c = 1.0 and down- 
stream of this location the isovels intersect the horizontal axis indicating 
the wake development. 
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Another view of the mean velocity field is shown in Fig. 5 which presents 

velocity profiles at various streamwise locations. The distances x/c and y /c 

s 

are defined as in Fig. 4 and u/u^^ is the local velocity component in the x- 
direction divided by the velocity at upstream infinity. The dimensionless velo- 
city reaches values greater than unity since the inviscid velocity at the edge of 
the boundary layer is greater than the velocity at upstream infinity. Although 
the flow approached a separated state at the trailing edge vicinity, separation 
never actually occurred; this could be due to the turbulent nature of the flow. 

The calculated turbulence kinetic energy contours are presented in Fig. 6. 

The maximum predicted value of turbulence kinetic energy, k = 0.06, occurs at 

x/c ~ 0.1 and very close to the airfoil surface where the assumed wall function 

for turbulence kinetic energy was applied. Since these high values occur both 

in the wall function region and in the region where the flow was transitional 

(Fig. 7) this maximum value should be viewed with some skepticism. Further 

downstream, the values of k in the vicinity of the wall decreased to values 

between 0.02 and 0.04 which are more in keeping with experiment (Ref. 61). 

Turbulent viscosity isobars are presented in Fig. 7. Figure 7 indicates that 

transition occurs in the region x/c 0.2 which is somewhat upstream of the 

region noted experimentally by Gregory and O'Reilly (Ref. 63). Further, 

according to Gregory and O’Reilly (Ref. 63) the transition location for flow 

£ 

about a NACA 0012 airfoil at Re = 10 is at x/c = 0.4. However, it must be 

c 

recalled that the present turbulence model was applied without any ’tuning’ 
and, hence the agreement is considered satisfactory for the present. As in 
Fig. 4, x/c = 1.0 represents the airfoil trailing edge and the wake is found 
downstream of this point. 

The second case considered is again for an airfoil at zero incidence, 
however, in this case flow about a full airfoil rather than flow about half an 
airfoil with use of symmetry conditions was calculated. A constructive grid 
having 81 x 30 grid points was used. The flow conditions for this second case 
were a Reynolds number based upon chord of 10^ and a Mach number of .147. In 
addition, the computational grid for this case was more highly stretched than 
for the previous case with the first grid point off the airfoil surface being 
located a distance of 0.0002 chords from the surface. When the full airfoil 
calculation is considered a major advantage of the constructive coordinate 
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grid becomes apparent. The Poisson generated coordinate system (Ref. 27) gives 
coordinate lines with discontinuous slope at the geometric centerline of the 
coordinate system emanating from the airfoil trailing edge, see Fig. 8, (the 
branch cut); the constructive coordinate system does not contain this problem. 
Therefore, performing a full airfoil calculation integration across the wake 
from the lower to upper boundary including the centerline point (e.g., line a-b 
on Fig. 8) presents no difficulty with the present coordinates. Although such 
a procedure is possible in a Poisson generated coordinate system, its implemen- 
tation would require using a device such as one-sided differencing at the 
centerline or not solving the equations on the centerline but determining values 
there by extrapolation. One-sided differencing may lead to numerical 
difficulties and the extrapolation technique is equivalent to treating the 
centerline as a boundary and applying boundary conditions along this interior 
line. With either of these two devices the centerline would be treated in a 
special manner and, therefore, both of these devices are expected to be inferior 
to simply treating the centerline points as interior points in the manner of 
the present analysis. 

The full airfoil calculation was made with function conditions applied for 
density and the two velocity components along the outer boundary line d-e-f-g-h. 
As for the symmetric half airfoil calculation the values were taken as the flow 
field values appropriate for a Joukowski airfoil of equal thickness. As in the 
first calculation the outer boundary was approximately three chords away from 
the airfoil surface. No slip conditions and zero density gradient conditions 
were applied at the airfoil surface and zero first derivative boundary con- 
ditions along the ' streamwise' direction for both velocity and density were 
applied at the outflow boundary. 

The results of the full airfoil calculation are very similar to those of 
the half airfoil calculation presented in Figs. 3-7. The flow field is 
predicted to be symmetric about the centerline as is expected. The surface 
static pressure distribution shown in Fig. 9 is only slightly different from 
that shown in Fig. 3. The major difference being the higher pressure in the 
suction surface trailing edge region; this has been discussed previously and 
may be a result of the better boundary layer resolution in the full airfoil 
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calculation. Likewise the velocity, turbulence energy and turbulent viscosity 
fields differed little from those previously presented. 

The third example considered is a NACA 0012 airfoil at 6° incidence. Again 
the Reynolds number based upon chord was 10^ and the Mach number was 0.147. 

As in the previous example an 81x30 grid was used. In running the calculation 
some difficulty was experienced in bringing the transverse velocity, w, 
smoothly into the outer specified function boundary condition in the region 
d-e ' and g f -h of Fig. 8. The initial attempt to resolve this problem focused 
upon the transverse grid which was very widely spaced in the vicinity of the 
outer loop, d-e’-f-g'-h. A new grid was constructed which sacrificed some 
minor amount of resolution in the boundary layer and even more resolution In 
the flow midregion to obtain better resolution near the outer loop. In this 
grid the first grid point off the airfoil was located a distance 0.00033 from 
the airfoil surface. Although this procedure decreased the oscillations in 
the w-component of velocity, it did not eliminate them; hence a second modifi- 
cation was made. In this modification the function boundary conditions on p, 
u and w along lines d-e and g-h (see Fig. 8) were replaced by specified static 
pressure and zero first derivatives of u and w; this set of boundary conditions 
removed the observed oscillations and the calculation proceeded without 
difficulty. Further this new set of boundary conditions appear more desirable 
from physical considerations. It should also be noted that the flexibility 
of the constructive coordinate in obtaining some desired mesh proved a very 
desirable property. 

The predicted pressure distribution is compared with the data of Gregory 
and O f Reilly (Ref. 63) taken for a Reynolds number 2.8x10^ in Fig. 10. As 
shown in Fig. 10 the major discrepancy between data and analytic prediction 
occurs in the leading edge region where the analysis fails to predict the 
correct suction peak. This discrepancy is at least partially a result of grid 
resolution. The strong favorable pressure gradient regipn leading to the 
suction peak occurs in a very limited region of the flow field between x/c 0 
and x/c *=» 0.01. This region extends over only one percent of the airfoil chord 
and only one tenth of one percent of the entire grid extent. In interest of 
computer run time economy the grid was limited to 81x30 points (a total of 2430 
grid points) and even though points were packed into the leading edge region. 
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only four pseudo-radial lines were placed within the favorable pressure grad- 
ient region. Thus even with a total of 2430 grid points and severe leading 
edge grid packing, resolution in this region was marginal. It is expected 
that increased resolution would result in better agreement with the data. 

A second area of moderate disagreement occurs in the trailing edge region; 
this results from differing predictions of trailing edge pressure. The 
discrepancy between predicted and measured trailing edge pressure also occurred 
in the zero incidence case and has been discussed previously in some detail. 

The lift distribution presented in the lower portion of Fig. 10 shows excellent 
agreement between data and analysis except in the region of the suction peak 
where as mentioned above the probable cause of disagreement is grid point 
resolution. 

In regard to other aspects of the flow field the predicted suction surface 

transition location occurs at x /c ?»0.08. The data of Gregory and O'Reilly 

gives x_ ps 0.04 for a Reynolds number of 2.8 x 10 and x /c «0.08 for a 

6 ^ 

Reynolds number of 1.48 x 10 . Thus the predicted transition location is in 
excellent agreement with data. The transition location predicted on the pressure 
surface is x^/c ~0.30; thus the pressure surface boundary layer has a consider- 
ably longer laminar run than does the suction surface boundary layer. This 
result is as expected. The flow on the airfoil suction surface is subject to 
a very strong acceleration as it proceeds from the stagnation point around the 
airfoil leading edge and the large acceleration with resultant generation of 
large velocities is shown in Fig. 11. As the flow proceeds downstream, it 
encounters an adverse pressure gradient leading to a thickening of the boundary 
layer, typical adverse pressure boundary layer profiles and eventually separation 
at x/c *=».85. The appearance of separation is at variance with the data of 
Gregory and O'Reilly (Ref. 63) taken at slightly higher Reynolds number which 
shows no separation at this incidence angle. However, the predicted separation 
region occurs on the very aft part of the airfoil and remains very small, thus 
having only a minor effect upon the generated pressure distribution. The 
difference between the analytic prediction of separation and the data may be due 
to any of several possibilities. Since the calculation was run at a slightly 
lower Reynolds number than the data, the difference may be due to Reynolds number 
effects. However, it is more likely that the discrepancy is tied to the 
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turbulence model; modification of the turbulence model structural functions 
could delay separation. Nevertheless, despite this difference the comparison 
between analytic prediction and experimental data is good. As is expected, the 
flow on the pressure surface is much less dramatic than that on the suction. 

The pressure surface has no suction peak and the boundary layer on the surface 
never encounters a strong adverse pressure gradient. The result is a well- 
behaved flow with no approach to separation. The suction surface and pressure 
surface boundary layers merge at the trailing edge where the wake is seen to be 
highly asymmetric. 

A detailed picture of the flow in the leading edge region is presented in 
Fig. 12. In this figure the splitting of the flaw as it encounters the airfoil, 
the appearance of the front stagnation and the large acceleration as the flow 
proceeds around the leading edge are all clearly pictured. 

The turbulence energy field on both the suction and pressure surfaces is 
shown in Fig. 13. Since the regions of high turbulence energy are concentrated 
very near the airfoil surface, they are presented as lines of constant k on a 
plot of distance along airfoil surface to distance from airfoil surface as was 
done in Figs. 4, 6 and 7. The field shows distinctly different characteristics 
on the two different surfaces. On the suction surface large turbulence energies 
are generated in the region of the high adverse pressure gradient. However as 
the boundary layer proceeded in the high adverse pressure gradient region and 
approached separation, x/c > 0.4, the turbulence energy drops as a result of 
the low transverse velocity gradients associated with nearly separated boundary 
layers which can no longer produce turbulent energy. In fact on the aft portion 
of the airfoil suction surface the maximum turbulence energy appears' away from 
the wall. The turbulence energy field on the pressure surface exhibits a 
different character. Although the field here does not generate as large k’s 
as on the suction surface, relatively high turbulence levels are present over 
the entire surface flow field. Again this is to be expected since the 
boundary layer in this region does not approach separation. 
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CONCLUDING REMARKS 


The present report describes a solution of the compressible, turbulent 
isolated airfoil flow field problem. The analysis used solves the full 
compressible Navier-Stokes equation via the consistently split linearized 
block implicit procedure of Briley and McDonald in a body-fitted constructive 
coordinate system. The constructive system has the airfoil surface as a 
coordinate line by definition and is capable of obtaining high grid resolution 
at user-specified regions of the flow field. To the author's knowledge the 
present results represent the first solution of the airfoil flow problem using 
the Navier-Stokes equations in conjunction with a constructive coordinate 
system. 

Several forms of the Navier-Stokes equations are presented as possible 
candidates to obtain an accurate and efficient solution. Based upon a variety 
of considerations, including sample calculations of low Reynolds number flow 
about a circular cylinder, it is determined that the strong conservation form 
of the equations should not^be solved; rather a form in which the geometric 
factors are not in conservative form is deemed superior for this problem. 

Since high Reynolds number flow is the flow of interest, the present effort 
proposes a turbulence model based upon the turbulence kinetic energy equation 
with structural functions specified so as to be valid for laminar, transitional 
and fully turbulent flow. When used in conjunction with the mean flow equation, 
the model predicts a turbulence kinetic energy field, as well as a mean flow 
field. The analysis has been applied to high Reynolds number, subsonic Mach 
number flow about a NACA 0012 airfoil. At zero incidence the analysis gives 
a good comparison between the predicted airfoil static pressure distribution 
and experimental data. The only region which contains a discrepancy is the 
trailing edge region where grid resolution or turbulence modelling effects may 
give larger viscous displacement effects than were observed experimentally. 

In addition, the mean flow and turbulence energy results are as expected. 

Regions of laminar flow appear in the airfoil leading edge region and in regions 
distant from the airfoil. Transitional flow occurs near the airfoil surface in 
the quarter chord region and turbulent flow is predicted over the aft part of 
the airfoil and in the airfoil wake. 
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At 6 deg. incidence agreement between predicted and measured pressure 
distributions was again good except in the leading edge region where grid 
resolution problems prevented an accurate prediction of the suction peak. 

The predicted suction surface transition point was in excellent agreement 
with data and the mean and turbulence fields were qualitatively as expected. 

Although some problems remain, particularly in regard to grid resolution 
and turbulence modelling, the results of the present analysis clearly show 
the capability of predicting the high Reynolds number isolated airfoil flow 
field in a constructive coordinate system including turbulence effects. 


35 


APPENDIX - THE NUMERICAL METHOD 


Linearization Technique 

A number of techniques have been used for implicit solution of the following 
first-order nonlinear scalar equation in one dependent variable £(x,t): 


=F(<£) 3g(^)/3x 


(Al) 

Special cases of Eq. (Al) include the conservation form if F(<J>) * 1, and quasi- 
linear flow if G(<f>) * <f>. Previous implicit methods for Eq . (Al) which employ 
nonlinear difference equations and also methods based on two-step predictor- 
corrector schemes are discussed by Ames (Ref. 64, p. 82) and von Rosenburg 
(Ref. 65), p. 56). One such method is to difference nonlinear terms directly 
at the implicit time level to obtain nonlinear implicit difference equations; 
these are then solved iteratively by a procedure such as Newton's method. 
Although otherwise attractive, there may be difficulty with convergence in the 
iterative solution of the nonlinear difference equations, and some efficiency 
is sacrificed by the need for iteration. An implicit predictor-corrector 
technique has been devised by Douglas and Jones (Ref. 66) which is applicable 
to the quasilinear case (G “ <P ) of Eq . (Al) . The first step of their procedure 
is to linearize the equation by evaluating the non-linear coefficient as F (d> n ) 
and to predict values of using either the backward difference or the 

n+i 

Crank-Nicolson scheme. Values for 4> are then computed in a similar manner 
using F(c|) n+ ^^) and the Crank-Nicolson scheme. Gourlay and Morris (Ref. 67) 
have also proposed implicit predictor-corrector techniques which can be applied 
to Eq. (Al) . In the conservative case (F*l) , their technique is to define 
G ( 4>) by the relation G(<J>) * 4>G ( tp) when such a definition exists, and to evaluate 
G(p n+ ^) using the values for computed by an explicit predictor scheme. 

With G thereby known at the implicit time level, the equation can be treated as 
linear and corrected values of 4> n+ ^ are computed by the Crank-Nicolson scheme. 

A technique is described here for deriving linear implicit difference 
approximations for nonlinear differential equations. The technique is based 
on an expansion of nonlinear implicit terms about the solution at the known time 
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level, t n , and leads to a one-step, two-level scheme which, being linear in 
unknown (implicit) quantities, can be solved efficiently without iteration. 
This idea was applied by Richtmyer and Morton (Ref. 68, p. 203) to a scalar 
nonlinear diffusion equation. Here, the technique is developed for problems 
governed by a nonlinear equations in z dependent variables which are 
functions of time and space coordinates. Although the present effort concen- 
trates upon two spatial dimensions and time, the technique will be described 
for the three-dimensional, unsteady equations. 

The solution domain is discretized by grid points having equal spacings 

12 3 12 3 

in the computational coordinates. Ay , Ay and Ay in the y , y and y 

directions, respectively, and an arbitrary time step. At. The subscripts 

12 3 

i, j, k and superscript n are grid point indices associated with y , y , y 
and t, respectively, and thus 4>? . , denotes 'Ky'f, y?, y?, t n ) . It is 

x , j , Ic l ^ j ic 

assumed that the solution is known at the n level, t , and is desired at 

the (n+1) level, t n+ ^. At the risk of an occasional ambiguity, one or more 

of the subscripts is frequently omitted, so that 4> n is equivalent to <J>? . . 

i » J » J 

Although present attention is focused on the compressible Navier-Stokes 
equations, the numerical method employed is quite general and is formally 
derived for systems of governing equations which have the following form: 


dH(4>)!d\ = 5b (<P) +S(<£) 


(A2) 

where 4> is a column vector containing dependent variables, H and S are 
column vector functions of , and 5b is a column vector whose elements are 
spatial differential operators which may be multidimensional. The generality 
of Eq. (A2) allows the method to be developed concisely and permits various 
extensions and modifications (e.g., noncartesian coordinate systems', turbulence 
models) to be made more or less routinely. It should be emphasized, however, 
that the Jacobian 3H/S<|» must usually be nonsingular if the ADI techniques as 
applied to Eq. (A2) are to be valid. A necessary condition is that each 
dependent variable appear in one or more of the governing equations as a time 
derivative. An exception would occur if for instance, a variable having no time 
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derivative also appeared in only one equation, so that this equation could be 
decoupled from the remaining equations and solved ji posteriori by an alternate 
method. As a consequence, the present method is not directly applicable to the 
incompressible Navier-Stokes equations except in one-dimension, where ADI 
techniques are unnecessary. For example, the velocity-pressure form of the 
incompressible equations has no time derivative of pressure, whereas the 
vorticity-stream-func tion form has no time derivative of stream function. For 
computing steady solutions, however, the addition of suitable "artificial" time 
derivatives to the incompressible equations, as was done in Chorin's (Ref. 69) 
artificial compressibility method, would permit the application of the present 
method. Alternatively, a low Mach number solution of the compressible equations 
can be computed. 

The linearized difference approximation is derived from the following 
implicit time-difference replacement of Eq. (A2) : 

(H n + , -H n )/AI= / 8[2>(^ n + , )+S n + ']+(l-/3)fS(0 n ) + s n ] 

(A3) 

where, for example, H n+ ^ = H(<f> n+ ^ ) . The form of and the spatial differencing 

are as yet unspecified. A parameter 3(0 < B < 1) has been introduced so as to 
permit a variable centering of the scheme in time. Equation (A3) produces a 
backward difference formulation for 6 = 1 and a Crank-Nicolson formulation for 
B = 1/2. 

The linearization is performed by a two-step process of expansion about the 
known time level t n and subsequent approximation of the quantity ( 3<$>/3 t) n At , 
which arises from chain rule differentiation, by (<t> n+ ^-<t> n ) . The result is 

H n+I = H n +(dH/d<£) n (<£ n + , -<£ n ) + 0 ( At ) 2 (A4a) 

S n + l = S n +Os/d^) n (</> n + , ~<£ n ) +0( At ) 2 (A4b) 
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2>W> n+l ) = ^(O+(^/d*K* n+, -* n ) + 0(At) z 

(A4c) 

The matrices 3H/3<}> and 3S/3<£ are standard Jacobians whose elements are 

defined, for example, by (3H/3<fi) = 3H /3<f> . The operator elements of the 

qr q r 

matrix 3^/3<J> are similarly ordered, i.e., (3.5>/3cJ>) H 3.$) /3<f> ; however, the 

qr q r 

intended meaning of the operator elements requires some clarification. For the 
q fc row, the operation (3^^/34>) n ( 4 > n+ ^ - <J> n ) is understood to mean that 
{3/3t^[4>(x,y,z,t) ] } n At is computed and that all occurrences of ( 34> r / 3 1 ) n 
arising from chain rule differentiation are replaced by - $^)/At. 

After linearization as in Eqs. (A4) , Eq. (A3) becomes the following linear 
implicit time-differenced scheme: 

ldH n /d<f>W ' + ' -*") /At = .2>(-*> n ) +S n + P (d^/df + d s n /3<W> + '-<#>") 

(A5) 


Although H is linearized to second order in Eq. (A4), the division by At in 
Eq. (A3) introduces an error term of order At. A technique for maintaining 
formal second-order accuracy in the presence of nonlinear time derivatives is 
discussed by McDonald and Briley (Ref. 70), however, a three-level scheme 
results. Second-order temporal accuracy can also be obtained (for g = 1/2) by 
a change in dependent variable to 0 £ H ( 4 > ) , provided this is convenient, since 
the nonlinear time derivative is then eliminated. The temporal accuracy is 
independent of the spatial accuracy. 


On examination, it can be seen that Eq. (A5) is linear in the quantity 
(4> n+ ^ - 4> n ) and that all other quantities are either known or evaluated at the 
n level. Computationally, it is convenient to solve Eq. (A5) for ($ n+ ^ - $ n ) 
rather than <p . This both simplifies Eq. (A5) and reduces roundoff errors, 
since it is presumably better to compute a small 0(At) change in an 0(1) 
quantity than the quantity itself. To simplify the notation, a new dependent 
variable ip defined by 


f = $-4> n 


(A6) 
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is introduced, and thus ■ <p n+ ^ - $ n , and =* 0. It is also convenient 

to rewrite Eq. (A5) in the following simplified form: 

lA + AtJ)*"* 1 = At [2> (*")+S n ] 


(A7a) 


where the following symbols have been introduced to simplify the notation: 


As dH n /d<t> ~/3&\{dS n /dcp) 


(A7b) 


2=l -Pid 2>/d<j>) 


(A7c) 


It is noted that is a linear transformation and thus H. (0) = 0. Further- 

more, if J2>(<£) is linear, then J! { ip) =* -ft Ib(ip) . 

Spatial differencing of Eq. (A7a) is accomplished simply by replacing 

i 2 i i 

derivative operators such as 3/3y , 3 /3y 3y by corresponding finite differ- 
ence operators, D^, D^. Henceforth, it is assumed that lb and J!. have been 
discretized in this manner, unless otherwise noted. 

Before proceeding, some general observations seem appropriate. The fore- 
going linearization technique assumes only Taylor expandability, an assumption 
already implicit in the use of a finite difference method. The governing 
equations and boundary conditions are addressed directly as a system of coupled 
nonlinear equations which collectively determine the solution. The approach 
thus seems more natural than that of making ad hoc linearization and decoupling 
approximations, as is often done in applying implicit schemes to coupled and/or 
nonlinear partial differential equations. With the present approach, it is not 
necessary to associate each governing equation and boundary condition with a 
particular dependent variable and then to identify various ''nonlinear 
coefficients" and "coupling terms" which must then be treated by lagging, 
predictor-corrector techniques, or iteration. The Taylor expansion procedure 
is analogous to that used in the generalized Newton-Raphson or quasi-linearization 
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methods for iterative solution of nonlinear systems by expansion about a known 
current guess at the solution (e.g.. Bellman & Kalaba, Ref. 71). However, the 
concept of expanding about the previous time level apparently had not been 
employed to produce a noniterative implicit time-dependent scheme for coupled 
equations, wherein nonlinear terms are approximated to a level of accuracy 
commensurate with that of the time differencing. The linearization technique 
also permits the implicit treatment of coupled nonlinear boundary conditions, 
such as stagnation pressure and enthalpy at subsonic inlet boundaries, and in 
practice, this latter feature was found to be crucial to the stability of the 
overall method (Ref. 40). 

Application of Alternating-Direction Techniques 

Solution of Eq. (A7a) is accomplished by application of an alternating- 
direction implicit (ADI) technique for parabolic-hyperbolic equations. The 
original ADI method was introduced by Peaceman and Rachford (Ref. 72) and 
Douglas (Ref. 73); however, the alternating-direction concept has since been 
expanded and generalized. A discussion of various alternating-direction 
techniques is given by Mitchell (Ref. 74) and Yanenko (Ref. 39). 

The present technique is simply an application of the very general 
procedure developed by Douglas and Gunn (Ref. 41) for generating ADI schemes 
as perturbations of fundamental implicit difference schemes such as the 
backward-difference or Crank-Nicolson schemes. 

For the present, it will be assumed that ,2)(<J>) contains derivatives of 

12 3 

first and second order with respect to y , y and y , but no mixed derivatives. 

In this case, lb can be split into three operators, associated 

12 3 i- Z J ■ 

with the y , y and y coordinates and each having the functional form 

2>, - 3/3y^, 3 ^/3y^3y^) . Equation (A7a) then becomes 

[a + At( +J 2 + )]* n + , = At[(3>, +^ 2 + ^ 3 )^> n + s n ] 

(A8) 

Recalling that J!( ip n ) = 0, the Douglas-Gunn representation of Eq. (A8) can be 
written as the following three-step solution procedure: 
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(A+At^ ty* = At[( +^ 2 + * 3 )*%S"] 

(A9a) 

(A + At J 2 )*** = At * 

(A9b) 

(A + At =A*** 

(A9c) 


where if) and 4> are intermediate solutions. It will be shown subsequently 
that each of Eqs. (A9) can be written in narrow block-banded matrix form and 

A AA 

solved by efficient block-elimination methods. If ij/ and ^ are eliminated, 
Eqs. (A9) become 

(A + At^.U-'fA + At y 2 )A"'(A + AI J )*" + ' = At [( 2 >, + 3> 2 + 2 > i ) ^ + S n ] 

(A10) 

If the multiplication on the left-hand side of Eq . (A10) is performed, it becomes 

2 

apparent that Eq. (A10) approximates Eq. (A8) to order (At) . Although the 
stability of Eqs. (A9) has not been established in circumstances sufficiently 
general to encompass the Navier-Stokes equations, it is often suggested (e.g., 
Richtmyer & Morton, Ref. 68, p. 215) that the scheme is stable and accurate under 
conditions more general than those for which rigorous proofs are available. 

This latter notion was adopted here as a working hypothesis supported by favor- 
able results obtained in actual computations (e.g.. Ref. 38). 

A major attraction of the Douglas-Gunn scheme is that the intermediate 

. * aa , , n +l 

solutions 4* and 4* are consistent approximations to V . Furthermore, 

n ^ ivH 

for steady solutions, 4< = 4* = 4* =4' independent of At. Thus, physical 

boundary conditions for 4> n+ ^ can be used in the intermediate steps without a 
serious loss in accuracy and with no loss for steady solutions. In thiB respect, 
the Douglas-Gunn scheme appears to have an advantage over locally one-dimensional 
(LOD) or "splitting" schemes, and other schemes whose intermediate steps do not 
satisfy the consistency condition. The lack of consistency in the intermediate 
steps complicates the treatment of boundary conditions and, according to Yanenko 
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(Ref. 39, p. 33), does not permit the use of asymptotically large time steps. 

It Is not clear that this advantage of the Douglas-Gunn scheme would always 
outweigh other benefits which might be derived from an alternative scheme. 
However, since the ADI scheme can be viewed as an approximate technique for 
solving the fundamental difference scheme, Eq. (A7a) , alternate techniques can 
readily be used within the present formulation. 

It is worth noting that the operator can be split into any number of 

components which need not be associated with a particular coordinate direction. 
As pointed out by Douglas and Gunn (Ref. 41) , the criterion for identifying 
sub-operators is that the associated matrices be "easily solved" (i.e., narrow- 
banded). Thus, mixed derivatives can be treated implicitly within the ADI 
framework, although this would increase the number of intermediate steps and 
thereby complicate the solution procedure. Finally, only minor changes are 
introduced if, in the foregoing development of the numerical method, H, .2), and 
S are functions of the spatial coordinates and time, as well as 4». 

Solution of the Implicit Difference Equations 

Since each of Eqs. (A9) is implicit in only one coordinate direction, the 
solution procedure can be discussed with reference to a one-dimensional problem. 
For simplicity, it is sufficient to consider Eq. (A9a) with = 0. 

Consider the following three-point difference formulas: 


O lr '#>s[aA_ + (l-a)A + ]*/A’Z = (a*/d*). + 0 IAX *+(0-1/2) Ax] 


(Alla) 


Oy^stA. A_)*/(ATt) Z = (d*+/d* 2 ) i + o(A‘3r*) 

(Allb) 

for a typical computational coordinate x. Here, A _ = ^ a + = 

and a parameter a has been introduced (0 < a < 1) so as to permit continuous 
variation from backward to forward differences. The standard central difference 
formula is recovered for a * 1/2. 
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As an example, suppose that the component of lb" has the form 


d 

+ F lV*>7~T ®2g( ^ ) 


(A12) 


where F and G are column vector functions having the same but an arbitrary 

T 

number of components; F denotes the transpose of F. The form of Eq. (A12) 
permits governing equations having any number of first and second derivative 
terms. Then, 


At 


at ~ a* 






aF, 

dj> 


r_ T d Z ^ G 2q ^ 2G 2q 

+ L 2Q aifz d<t> ~&r 







(A13) 


It is now possible to describe the solution procedure for Eq. (A9a) for 

the one-dimensional case with given by Eq. (A12) and appropriate difference 

x 2 

formulas. Because of the spatial difference operators, D~ and D~ , Eq. (A9a) 

* * * xx 

contains an< * consequently, the system of linear equations 

generated by writing Eq. (A9a) at successive grid points x^ can be written in 
block- tr id lagonal form (simple tridiagonal for scalar equations, l * 1) . The 
block-tridiagonal matrix structure emerges from rewriting Eq. (A9a) as 




+ b!V* 

ii *i+i 



(A14) 


where a, b, c are square matrices and d is a column vector, each containing 

only n-level quantities. _ When applied at successive grid points, Eq. (A14) 

* 

generates a block-tridiagonal system of equations for which, after appro- 

priate treatment of boundary conditions, can be solved efficiently using 
standard block-elimination methods as discussed by Isaacson and Keller (Ref. 75, 
p. 58). The solution procedure for Eqs. (A9b&c) is analogous to that Just 
described for Eq. (A9a) . It is worth noting that the spatial difference para- 
meter a can be varied with i or even term by term. For example, an "upwind 
difference" formula can be obtained if a is chosen as 1 or -1 depending on 
the sign of the elements of F^. 
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Figure 3. - Surface pressure distribution for NACA 0012 airfoil at zero incidence 
(half airfoil calculation). 
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Sketch of Poisson generated coordinate system 



Figure 8. - Comparison of coordinate system generated by 
solution of Poisson equation and constructive 
coordinate system. 
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Figure 9- - Surface pressure distribution for NACA 0012 airfoil at zero incidence 
(full airfoil calculation). 
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Figure 11.- Velocity field for NACA 0012 airfoil at 6° incidence, 
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